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1. Introduction 



^ . The numerical iiydro dynamical modelling of radial stellar pulsations 

goes back to the pioneering efforts of Christy (1966) in the mid 1960s. 
Reviews of the stellar pulsation problem can be found in the literature 
{e.g. Cox & Giuh 1969, Cox 1980, Buchler 1990). In our own past 
numerical work we have adopted the hydrodynamics code developed 
by Stellingwerf (1974, 1975) on the basis of Praley's (1968) implicit 
\ Lagrangean scheme. This is a very stable and robust code, and in its 

relaxation version it is an excellent tool for computing periodic pulsa- 
tions and their linear stability (Floquet analysis). The use of the code 
\ has produced good results for a variety of stellar models ranging from 

0^ ■ RR Lyrae {e.g. Stellingwerf 1975, Kovacs & Buchler Kovacs 1988a) and 

r~| ! Cepheids {e.g. Moskalik, Buchler & Marom 1992) to irregularly pulsat- 

0^' ing variables (Kovacs &: Buchler 1988b), and good overall agreement 

Q . with observations can be achieved. As these theoretical studies have 

|H \ progressed, and as more and better observational has accumulated it 

^ ■ has become increasingly clear that a better code is needed {e.g. Buchler 

1990) if further advances are to be made (Kovacs 1990, Kovacs & 
Buchler 1993). 

^ ■ In fact, right from the beginning of Cepheid modelling misgivings 

^ \ were expressed about the spatial resolution that can be obtained with 

a Lagrangean grid. The earliest adaptive attempts made to abate this 
deficiency were those of Castor, Davis &: Davison (1977) and Aikawa & 
Simon (1983), Simon & Aikawa (1986), but these codes were not energy 
conserving. More recently, taking advantage of progress in numerical 
mathematics, a couple of efforts specifically aimed at stellar pulsations, 
and parallel to ours have been made (Dorfi &: Drury 1987, Dorfi & 
Feuchtinger 1991, Feuchtinger & Dorfi 1994, Cox, Deupree & Gehmeyr 
1992, Gehmeyr 1992a,b, 1993). The three codes differ in many respects. 
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however. In particular, the other two codes solve the internal energy 
equation whereas we solve the total energy equation (While analytically 
these equations are equivalent, the tests that we describe below show 
that, numerically, the solution of the total energy equation is more 
accurate. Furthermore our grid equation while similar to Dorfi's uses a 
mass density rather than a radius density. A comparison of the three 
codes on a few standard stellar models would be desirable. 

Our strategy in developing an adaptive code has been very conserva- 
tive. Because of the success of the Stellingwerf-Praley scheme we have 
decided to build an adaptive code on top of this Lagrangcan scheme. 
By putting the mass and total energy equations in conservative form, 
and by properly defining the fluxes, we can write the adaptive scheme 
so that, globally, mass and total energy conservation are preserved. 
These constraints are essential for the stellar pulsation problem (c/. 
Buchler 1990; see also the tests below). Indeed, an examination of the 
energetics shows the pulsation to be rather delicate. The pulsation 
kinetic energy (~ 8 x 10^^ erg) is generally very small compared to 
both the gravitational energy (ss 8 x 10^^ erg) and the internal energy 
(!=a 4 X 10^^ erg) which have opposite signs (from the virial theorem 2K 
+ G ~ 0). In addition, in the nonadiabatic case, the growth-rates are 
often small, i.e. it takes many pulsations to build up to the saturation 
amplitude. In contrast, the heat that flows through the stellar envelope 
in one pulsation (si 2 x 10^^ erg) is very large (The indicated numerical 
values arc for a typical classical Cepheid) . The coupling of the pulsation 
to the heat-flux thus plays a very important role in the pulsations. It 
determines both the linear vibrational stability or instability and the 
saturation amplitude. Needless to say, it therefore requires an accurate 
difference scheme. One of the important checks is that in the linear 
regime the growth-rates obtained from the numerical hydrodynamics 
agree with those obtained from a linear stability of the normal modes 
of oscillation. 

The prescription for converting to an adaptive mesh is well known 
{e.g. Winkler, Norman 8z Mihalas 1984). Three velocities appear, namely 
u, the physical velocity of the fluid with respect to a spatially fixed 
reference frame (here, the equilibrium stellar model), Ugj-id, the veloc- 
ity of the mesh points with respect to the same reference frame and 
Urei = u — Ugrid the fluid velocity relative to the mesh. 

The adaptive mesh hydrodynamic equations, in their usual notation, 
are given by 




(1) 
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+ A [{pu + Fl) A]=0 



where 6 denotes the temporal change of a variable or expression and 
A the spatial variation. We treat the radiative flux in the diffusion 
approximation as we have done in the past with the Lagrangcan code. 

These equations have to be completed with a mesh equation which 
specifies how the mesh points are linked to the fluid. There is of course 
no universal prescription for a mesh equation, and we have adopted a 
variant of Dorfi & Drury (1987) which we shall describe in detail below. 

In a separate paper we will relax the equilibrium diffusion approxi- 
mation and add separate radiation equations. 



2. The Code 

Fraley's scheme is a staggered Lagrangean scheme in which the basic 
radii, velocities are defined at the cell edges (i), and the density and 
thermodynamic variables are defined at the cell centers as shown 

in Figure 1. 

Following common notation, each variable is assigned a spa- 
tial subscript i (or i + ^) and a temporal superscript n. Intermediate 
temporal values of the independent variables are defined with weights 
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Figure 1. The numerical scheme 



= + (1 - 0Jy" (7) 

The superscript (n + i ) on a function of independent variables, 
e.g. the area ^""'"2 , thus denotes the temporal average of the function. 
Spatial differences are defined by Ay, = Ui+i — Vi-i and ^Vi+i = 
Ui+i - Vi, respectively. 

Our 3''"^ equation is for the total energy {i.e. the sum of the internal, 
gravitational and kinetic energies) rather than for the internal energy 
alone as in Praley's scheme. However, Praley's scheme is such as to 
actually conserve the global total energy exactly. Therefore our differ- 
encing reduces to Fraley's, to within numerical round-off errors when 
the grid motion is Lagrangean {urei = 0). 

2.1. Continuity equation 
With the usual definitions 



the continuity equation takes on the form 

i^M^^+i - i^Mf^, + - Fl^} = 0, (8) 



where Fj^-^ is the mass flux through the zone boundary during the 
timestep St. 

There are many different ways of specifying the flux. We have chosen 
the second order scheme of van Leer (1979) for our code. (We have not 
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experimented with higher schemes, e.g. PPM; such schemes are much 
more cumbersome and expensive to implement, but arc not expected 
to increase the accuracy in the shocks and ionization fronts which 
are already well resolved, although they might increase the accuracy 
elsewhere) . In the following, the donor cell scheme is used only for the 
purposes of comparison. 

5t (pH.e/^)r^ = FmT' = ^t^ThpThu^^ (9) 

<^KeM^ = - (^r' - (10) 

with the switch S and dp are defined as 

S^{u) = l{l±sgn{u)) (12) 
dft+i = ^^^±i^5+(Ap,+iAp,). (13) 

For the purpose of comparison only we also use the less accurate first 
order upwind scheme (donor cell): 

{p"^^^) = ^~ {Urel,i)p^^i + 5+ iUrel,i)p^^i (14) 



2.2. Momentum Equation 
Again with the usual definition 

DMi = l(^DM,^^+DM,_^) (15) 
the two momentum source terms are differenced as 

Mp = 5t {AAp)i = 5t{Ar^i (p^^^ - p^+l) (16) 

^G = St(^^pAV^ =dtG{^)M^^^DM^^'^ (17) 

In the Fraley scheme global energy conservation is possible through 
the introduction of special averages. We therefore introduce the same 
averages in order to preserve this conservation in the Lagrangean limit 



3 

1 



(19) 
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For the fluxes we again use van Leer's expressions 
with 

Au,-, 1 Ait,_ 1 
dui = S+ (Aui, 1 An. 1 ) 

6Mi+^=i{SMi+^+6Mi). 
Finally the momentum equation becomes 

u^+^DM^+' - u'iDMl' + (^f;]+\ - +Mp + Mg = ^ (24) 

There are different ways of introducing 5Mi. Algebraically it is equal 
to Fj^^ , but numerically this form can cause problems (c/. comments 
by Winkler & Norman 1986, §V.B. on the use of 6M). The direct 
replacement of 5Mj by the mass flux also widens the band matrix in 
the iteration for each timestep as Eqs. 23 and 29 show. Both problems 
are avoided with the introduction of Mj as an independent variable and 
its deflnition 

AM,+ i =ft+iAV:,+ i (25) 
(^M"+^ = M^+^ - (26) 



(20) 

(21) 
(22) 

(23) 



2.3. Total energy equation 

A small complication arises because in Fraley's staggered scheme the 
internal energies are most naturally defined in the cell centers, but the 
kinetic and gravitational energies are defined at the cell boundaries. 
We therefore define the cell's total energy as being its internal energy 
plus one half of the combined kinetic and gravitational energies of its 
edges. 

^i+i = e^+i i^M^+i + gi+iDMi+i + giDMi (27) 
In accordance with Fraley's scheme we define the 'pressure-work' 
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{jmA)^'- = «+^) n^^ I (p:;^ + Pt?) (29) 
The total energy equation is 



(30) 

Note that the equation is written in fuh conservation form as a set of 
spatial differences, A. We define the relative energy flux again in the 
van Leer's fashion. Care has to be taken to center the corresponding 
fluxes correctly i.e. consistently with Eq. 27. First we calculate the 
advection terms of the kinetic energy at the cell centers, then calculate 
the spatial average at the cell boundaries where the internal energy flux 
is naturally defined. Other kinds of averaging, such as first calculating 
the the total energy density at cell centers and then advecting this 

quantity by SM^ ^ , can lead to numerical instabilities. 



FeP = St 



n+ 



n 

i+2 



1 



e, 

^9i+ 



->S+(Ae,+iAei 



i_ 3 — e,-_ 1 
1 Ao,_ 1 



(31) 

(32) 

(33) 

(34) 

(35) 



gi+i — gi-i 

The luminosity L is taken here in the radiation equilibrium diffusion 
approximation 



1 d aT^ 



(36) 



k{p, T) dm 3 

We generally use the difference form suggested by Stellingwerf (1974) 
which averages T^/n, but keep the option of using Stobie's (1969) 
average of the opacity. 

Finally a word about time-centering. Fraley's scheme is stable for 
> I . It is therefore tempting to increase the stability of the code by 
making all 6 > 0.5. However, in order to globally ensure total energy 
conservation in the Lagrangean limit, it is necessary to use 9u = ^ for 
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the velocity (Fralcy 1968). Furthermore, decentering the pressure in 
Eqs. (16) and (29) leads to improper mode growth or decay. In fact, 
it can be shown that even in the adiabatic limit the linear eigenval- 
ues incorrectly become complex (Kovacs, private communication); it is 
therefore important to set Op = ^ as well. Following Stellingwerf and our 
past experience we use ^ = | for the luminous flux L which has a small 
effect on the global results, except that is numerically stabilizing and 
reduces jitter. (Aside from the special Fraley averages {A) and (1/-R^) 
the luminosity L is therefore the only function that is not centered in 
time) . 

In some calculations, for test purposes only, we replace the the total 
energy difference equations (30) by the appropriate internal energy 
difference equations, viz. 

(t^r+t-t^r+i) + St[{L:^,^-L:^i)+p,^.A{{A)u),^,) (37) 

(38) 

where the internal energy is Ui_^_i = Cj+i AMj^i, the internal energy 

flux F^p = Stie'^^^dMp^ with(e"+^) as defined in Eq. 32. 

We note that this internal energy equation is also the equation solved 
by Gehmeyr and by Dorfi et al. {loc. cit). We note in anticipation that 
the performance of the code greatly deteriorates when this substitution 
is made. 

2.4. PsEUDO- Viscosity 

Finally a word about pseudo-viscosity. Our pressure everywhere in- 
cludes a von Neumann-Richtmyer (vNR) viscosity with a StelUngwerf 
cut-off parameter a^, involving the local sound velocity c, of the form 

Pvisc,i+^ = CqPi+1 ( \{Ui+i - Ui) - aqCg^i+i] ) (39) 



In the case where spherical geometric effects are important this expres- 
sion is advantageously replaced by the form suggested by Whalen (c/. 
Buchler & Whalen 1990, Eqs. 9-10). 

It may seem strange to still use pseudo-viscosity when the applied 
mathematical literature abounds in 'better' ways of treating shocks 
{e.g. Buchler 1990). First, we note that the use of flux limiters and 
Riemann solvers {e.g. Roe in Buchler 1990) assumes that the hydro- 
dynamics equations are written in a very specific 'conservative' form 
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which is not the case in our scheme, because of our desire to preserve the 
good features of the Fraley scheme. Eqs. 1-6 show that in a Lagrangean 
or near Lagrangean situation all but the 'pressure-work' flux and the 
radiative fluxes disappear and that the code would be unstable in the 
absence of pseudo-viscosity. Since, during the pulsation, parts of the 
mesh can indeed behave in a manner close to Lagrangean the possibility 
of such instability is not desirable. In our 'defense' we also note that 
when the grid motion redistributes mesh points into the shock, the 
effects of the pseudo- viscosity are largely reduced. 

2.5. The Mesh Equation 

Within constraints of stability, and ultimately accuracy, the choice 
of the mesh motion is to a large extent an art. On one extreme are 
the Eulerian and Lagrangean descriptions with a 'rigid' mesh, with 
Ugrid =constant and u^ei = 0, respectively, and on the other the mesh 
motions which follow the dynamics so rapidly and faithfully that the 
physical quantities are largely constant. A good mesh- function clearly 
lies in between. 

We have had good experience with Lagrangean meshes in which the 
outer zones, say up to hydrogen ionization (« 11,000 K) were equally 
spaced and the mesh size then increased geometrically inward. In the 
spirit of our conservative approach we therefore want to stay as close as 
possible to this type of mesh both in the static model and in the hydro- 
code. We want to use the adaptive features primarily, first, to resolve 
the partial hydrogen ionization during the pulsation, and second, to 
resolve shocks better than is possible with a standard Lagrangean mesh. 

For these reasons we adopt a mesh density defined in terms of the 
mass, rather than the radius, viz. 

1 = — (40) 

but for tests we keep the option of using the radius in the mesh density 
(i.e. nj_,_i = Xj+i /Ai?j_,_i ). The quantity X^_^_i is a mass scale (or 
length scale) . In the case of periodic pulsators our methodology will be 
to use first the fast Lagrangean StcUingwcrf code to obtain the limit 
cycles and then to switch over to an adaptive mesh with a minimum 
of transients. For this reason we use the cell masses AMf°l of the 

equilibrium model for the scaling, viz. Xi^i = Aikr^\ . This scale X^+i 

is thus held constant during the calculation. With this definition of 
1 the grid density is uniformly equal to 1 in the equilibrium model 
and m the Lagrangean code. We note that in general the magnitude of 
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depends on the definition of X^_^_i and really has a meaning only 
relative to the values of the grid density at different zones; in other 
words, generally the value of the grid density has nothing to do with 
the grid resolution. 

Otherwise we adopt the mesh diffusion (parameter a) and relaxation 
time (r) of Dorfi & Drury (1987): 

fii+i = n^+i - a(l - a)(nj+3 - 2^^+! + n^.i ) (41) 

For the structure function which gives the necessary grid resolution we 
have found it convenient to use the following function: 



(43) 



2 



where the scale factor 



Si,i=- /'"^ (44) 

l + (5n,+ i)4 ^ > 

is introduced to give an upper limit to the grid density, viz. 1/s. With- 
out the fourth order term in the denominator the grid density could 
increase unlimitedly causing numerical problems at velocity disconti- 
nuities: The artificial viscosity (Eq. 38) puts a constant number of grid 
points in the shock region, a number that is governed by the parameter 
Cq. This leads to constant velocity difference Ait independently of 

(k) / (k) 

the current grid density. If then the Af.'i /F.'i term in Eq. 42 is 

dominated by Au, the larger value of grid density generates a larger 
value of -^i+i and causes another increase of n. This positive feedback 
is stopped by the upper limit of grid density. 

Generally we introduce two terms in the sum (42), viz. 

Afll\=T,^s-T,_., <\=r,+ . (45) 

A/2\ = Ui+i - Ui, F^l\^ = Uo (46) 

Here the quantity Uq represents the velocity scale of the problem which 
can be mesh and time dependent. The /^^^ function serves the purpose 
of placing meshpoints in the regions of rapidly varying temperature, 
such as in the partial hydrogen ionization region which has a very 
sharp temperature variation, particularly in classical Cepheids. The 
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/(^) function acts only when strong velocity gradients are present and 

its obvious purpose is to resolve shocks. 

With these preambles the mesh equation finally is 

i^ = i^ (47) 

The specification of our mesh function is incomplete without bound- 
ary conditions. First we impose that the mesh be Lagrangean in the 
first few inner zones and in the last few zones. Second, following Dorfi h 
Drury (1987) we impose that nj_i = nj_|_i at both Lagrange-adaptive 
interfaces. 

With our definition of the grid density = AMf^\ /AMj_|_i ) the 

grid motion is Lagrangean when /J^ = because -^j+i = 1 and AMj_j_i 
remains constant. In other words the mesh equation is always satisfied 
when the (3k parameters are set to zero and the mass distribution is 
given by the initial Lagrangean model (-^j+i = '^j+i = !)• By a gradual 
increase of the (3k it is possible to evolve the Lagrange grid into the 
adaptive one (c/. also §3.3). 

For test purposes we have also experimented with the radius differ- 
ences in the grid density with = Ai?*"i (i.e.nj_,_i = Ai?*" i/Ai?j_,_i). 

In that case the (3k = condition gives a simplified grid equation: 
A-Rj^i /Ai?j_i =constant, which defines a quasi Eulcrian grid motion. 
Since several zones at the top and the bottom of the model are La- 
grangean, the grid motion of these special zones are given by the fluid 
motion: Rgr l — -^l- With these boundary conditions the solution of 
the simplified grid equation becomes 

where Rl^ and R^^ are the radius variations of the two Lagrangean 
zones bracketing the quasi Eulerian grid. 

3. Tests 

As performance tests for the code we have chosen the challenging Noh 
spherical shock problem (Noh 1987) and Sedov's (1959) point explo- 
sion problem. Analytical solutions exist in both cases. A perfect gas 
equation of state with (7 = 5/3) has been used for these tests. We will 
also present calculations for Cepheid model pulsations with a realistic 
equation of state and opacities to show the application of the code to 
stellar oscillations. 
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Figure 2. Numerical solution of the spherical Noh problem at t=0.6. Upper panels: 
tensor viscosity, lower panels: vNR viscosity, AM is used in the mesh density. 
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Figure 3. Spherical Noh problem. AR is used in the grid density; vNR viscosity. 
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3.1. Spherical Noh Problem 



The initial conditions for the Noh problem are uniform velocity {uj=- 
1) and density {pj=l) with zero pressure and temperature. The shock 
generated at the center of the sphere moves outward with a constant 
speed (t;=l/3). The density has a fiat profile behind the shock with 
p=64 at all times. 
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Figure 4- Spherical Noh problem. The internal energy eq. is solved instead of the 
total energy eq.; vNR viscosity 



For the numerical solution we have used a mesh with 50 zones 
initially spaced with equal AR. Since the spherical geometric effects 
are important we have also used the tensor artificial viscosity for these 
calculations (For a comparison of different treatments of artificial vis- 
cosity in Lagrangean calculations see Buchler and Whalen 1990). The 
parameters of the calculations are the following: (7^=1.0, = 1000, 
5=0.0002, r=0.0005, a=2.0 and Uo = 1. 

Figure 2. displays the results of the calculations for different pa- 
rameters together with the analytic solution. For the upper and lower 
row the Whalen tensor and the von Neumann-Richtmycr viscosity have 
been used, respectively. The adaptive calculation with the first order 
upwind scheme (second column) considerably improves the solution 
over the Lagrangean limit (/?„=0) which exhibits unphysical oscillations 
behind the shock. The use of the van Leer advection provides superior 
results (third column). This test clearly demonstrates that the higher 
order advection scheme is preferable in the adaptive calculation. The 
results arc independent of the choice of the artificial viscosity, in fact the 
use of the tensor artificial viscosity is not necessary with the adaptive 
computation even in this problem with strong sphericity. 

In Fig. 3. we present the solutions with an alternate form for the 
grid density, viz. with AR^^i in the denominator in Eq. 39. The vNR 
artificial viscosity is used. The first panel shows the quasi Eulerian 
result {f5u=0). The result with AR^^i is a little better than the one 
with AMj_,_i , especially at the center, l)ut the difference is insignificant 
elsewhere. (The glitch near the center is due to the inaccuracies intro- 
duced by the steeply varying cell mass in the initial profile which has 
been constructed with equal AR.) 

We have to emphasize that solving the total energy equation (Eq. 30) 
has been crucial in getting a correct solution for the Noh problem. We 
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Figure 5. Numerical solution of the Sedov problem. Prom left to right: initial velocity 
profile, Lagrangean and adaptive results. 



illustrate this in Figure 4 which shows the results for a calculation 
in which we now solve the internal energy equation (Eq. 437). It is 
true that with the adaptive mesh code we still get a sharp shock, but 
the velocity of the shock front (and thus the density) is very different 
from the correct values. When we monitor the time-dependence of the 
energies we notice a small gradual decrease in the global total energy 
(4.2% at t=0.6). Furthermore, the kinetic energy is too large and cor- 
respondingly the internal too small by ~ 7.3% and 15.6% respectively, 
whereas with the solution of the total energy equation the discrepancy 
in kinetic energy is less than 0.5%, i.e. 10 times less. In the Lagrangean 
limit, of course, the two energy equations yield the same solution to 
machine accuracy (because of the already mentioned property of the 
Fraley scheme). 

The order of the advection scheme plays a much more important role 
in the accuracy of the calculations when the internal energy equation 
is solved, as a comparison of Fig. 2 (top middle and top right) and 
Fig. 4 (left and middle) shows. This is an indication that the advection 
errors are much larger for the internal energy equation. An additional 
reason for the inaccuracies is found in the differencing of the source 
term (Eq. 37) pA{{A)u) which can vary very rapidly throughout the 
shock. In contrast, the total energy equation contains only flux terms 
and the corresponding Hugoniot-Rankine equation is automatically 
satisfied through near-discontinuities. 

The middle and right graphs of Figure 4 show that the inaccuracies 
associated with the solution of the internal energy equation occur for 
both definitions of mesh density (mass vs. radius). 
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3.2. Sedov Problem 

The analytic solution of an intense point explosion in a power-law 
density gradient with spherical symmetry is given by Sedov (1959). The 
motion of the gas is self-similar i.e. the shapes of the pressure, density 
and velocity profiles remain the same during the explosion. For an ideal 
gas with 7 = 5/3 the distance of the shock front from the center (Rs) is 
proportional to t^/^, and the physical parameters behind the shock are 
given by the initial value of the pressure, density and shock velocity. 
For an initial density given by p{R) = AR~'^ the solution has the form 

R 4AR A / rY 

for radii R less than the radial distance of the shock (Rg). The energy 
of the explosion appears in these functions only through the location 
of the shock front (Rs) for a given time. 

We choose the explosion energy so that the position of the shock at 
t=l is Rs=5 (in arbitrary units). Our initial model is constructed by 
discretizing the exact solution for t=l. Half of the 50 zones are placed 
below the front. The zones arc initially equally spaced in AR, but with 
different spacings below and above the front to get approximately the 
same zone mass on both sides. We use AM in the grid density. The 
numerical parameters are Cq=1.0, = 1000, s=0.005, r=0.003, a=2.0 
and Uo = I. 

The Lagrangean and adaptive solutions at t=8 are shown in Figure 
5 for the van Leer scheme. For comparison only Donor Cell results are 
also shown. By this time the intense shock, with density ratio 4:1 across 
the front, has progressed to a density 1/16 times the density at the start 
of the calculation. 

When the initial grid distribution does not satisfy the mesh equa- 
tion, which is the general situation, transient oscillations appear in 
the solution. Figure 5 shows that the initial transient leaves a large 
incorrect permanent scar on the solution when the inaccurate donor cell 
advection scheme is used. Only when we use the second order scheme 
are the fluctuations in the velocity minimal, and the position of the 
shock front is correctly obtained. 

As for the Noh problem, we again find it important to solve the total 
energy equation rather than the internal energy equation. With the use 
of the internal energy equation, the resulting shock is still sharp, but it 
lags. For example, at t=8 the position of the front is i?s=19.4 instead of 
the correct value of Rs=20. The accumulated error in the total energy 
reaches 10% for this time. The total kinetic energy which should be 
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constant during the shock propagation is also decreasing at the same 
rate. It is this underestimation of the kinetic energy that causes the 
error in the shock position. 

3.3. Cepheid Pulsations 

Table I. Typical errors of the total energy - Cepheid 
model problem. 





TUl) [K] 


^A{Etot) 


e(Etot)/P 


Scheme 


60 


5 10* 


4.0 10~^ 


-9 


-1.0 10"^ 


van Leer 


40 


1 10^ 


1.5 10" 


■7 


-1.5 10"® 


van Leer 


30 


2 lO'^ 


4.0 10" 


■7 


+1.0 10"* 


van Leer 


60 


5 10* 


3.0 10^ 


-8 


-6.0 10~* 


Donor Cell 


40 


1 10^ 


4.0 10^ 


-6 


-5.0 10-'^ 


Donor Cell 



We have selected a model from one of the Cepheid sequences of 
Moskalik, Buchler and Marom (1992). We use the OPAL opacities 
(Iglesias et al. 1992) matched with the opacities of Alexander and 
Ferguson (1994) at low temperatures. The chemical composition is 
X=0.7, Z=0.02. The computations are performed with 120 mass zones. 
The static model has been built as follows: equal mass zones inward 
from the surface; a temperature anchor of 11,000K in the 30th inward 
zone; the mass of the inner zones increases geometrically up to an inner 
temperature of Tj„ = 2 x 10^ K (This is the same type of grid that has 
been found to give good results in the Lagrangean code, e.g. Moskalik 
et al. 1992). This mass distribution also gives the reference grid for the 
adaptive calculations by setting the values of X^_^ i , as discussed earlier. 
With the stellar parameters M=6.5Mq, L=7213L0 and re//=5404 K 
the period of the model is ^ 18 days. 

It is notoriously difficult to get rid of long lasting transient oscilla- 
tions in an adaptive code. We have found the following procedure to be 
very useful in the case of periodic pulsations (asymptotic limit cycle) . 
A limit cycle solution is first calculated with the relaxation method of 
Stellingwerf applied to the Lagrangean code (starting from a perturbed 
static equilibrium model). This pulsating model is then inserted into 
the adaptive code in its Lagrangean mode, and the adaptive mesh 
is gradually switched on. More specifically, the (3^ parameters of the 
structure function are increased during the first s periods according to 
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the following recipe: 

Pk = 10^*/^^-=^/?^, (50) 

and they are then held constant for t > sP. This smooth switching 
on of the grid parameters transforms the Lagrangean mesh into an 
adaptive one. Since the grid equation is satisfied at the beginning of 
the calculation {i.e. in the /3fc=0, Lagrangean limit) we can evolve the 
model from Lagrangean to adaptive smoothly by this method and avoid 
annoying and often very slowly decaying transients. When the advec- 
tion errors are small (see below) the code converges from a Lagrangean 
to an adaptive grid limit cycle in several tens of periods. 

The transition of the grid from Lagrangean to adaptive is shown 
in Figure 6 for s = 2. The positions of the zones are given by M — 
Mj (j=70-120). The location of the partial ionization region is clearly 
revealed by the higher grid density which it attracts during the second 
period of pulsation. 
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Phase Phase 

Figure 7. Temperature profiles for the bump-Cepheid model Left box: Lagrangean 
grid, right box: adaptive grid 



The temperature parameter of the structure function is set to Pj^ =100 
For most of the tests we have used (3° =0 (the resolution in temperature 
is enough to get a smooth variation) , but we also have made calculations 
with (3° =10 and Uq = lO^cm/s. This results in a sharper shocks, 
but too large a value of P° attracts too many points into the shock 
and away from the partial ionization front, so that the latter loses its 
sharpness. The given values represent a good compromise between the 
two tendencies for this number (120) of total mesh points. For the other 
parameters we have chosen the values s=0.005, r=0.1 days and a=2.0. 
The pseudo viscosity in the form of Eq. 43 is used with Cq=4 (for some 
tests Cq=l) and 0^=0.01. 

The spatial switch-over from Lagrangean to adaptive mesh occurs 
at zone Jl ■ When it is made too far inside the stellar model (small Jl) 
we experience a numerical instability. This happens whether we solve 
the total energy equation or the internal energy equation. We note that 
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Phase Phase 

Figure 8. Velocity profiles for the bump-Cepheid model. Left: Langrangean, right: 
adaptive grid. 




Figure 9. Light-curve. Left: Langrangean, right: adaptive grid. 
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Table II. Fourier parameters of the limit cycle solution - radius variation. 
(LA: Lagrangean; TE: total energy, IE internal energy Eqs. are solved; vL: 
van Leer, DC: donor cell. 









Scheme 


Po 


Al 


R21 


R31 


ct>21 


031 


- 


4 


LA 


- 


17.6 


22.6 


0.193 


0.017 


6.159 


0.725 


- 


1 


LA 




IT.G 


2:5.2 


{).1!)7 


{).()2.1 


(i.lTG 


0.906 


60 


4 


TE 


vL 


17.6 


22.5 


0.197 


0.026 


6.190 


0.534 


60 


4 


IE 


vL 


17.6 


22.5 


0.195 


0.025 


6.194 


0.569 


60 


1 


TE 


vL 


17.6 


23.1 


0.199 


0.027 


6.200 


0.824 


60 


1 


IE 


vL 


17.6 


23.1 


0.199 


0.027 


6.203 


0.842 


40 


4 


TE 


vL 


17.5 


22.1 


0.186 


0.024 


6.213 


0.590 


40 


4 


IE 


vL 


17.5 


21.3 


0.182 


0.017 


6.177 


0.307 


40 


4 


TE 


DC 


18.7 


35.4 


0.347 


0.159 


6.265 


6.121 


40 


4 


IE 


DC 


18.0 


17.0 


0.148 


0.070 


0.395 


1.176 


30 


4 


TE 


vL 


17.3 


17.5 


0.138 


0.033 


0.127 


1.207 


30 


4 


IE 


vL 


17.3 


16.6 


0.130 


0.018 


0.021 


1.185 



the same problem has been encountered by Fcuchtingcr &: Dorfi (1994). 
The reason ultimately lies in the mesh equation. The global nature of 
the mesh equation and the mesh motions that occur in the outer zones 
because of moving physical features (mostly the H ionization front) 
induce small mesh motions even in the relatively quiescent inner zones. 
However, because of the large masses of these zones the internal and 
total energy advection errors can therefore reach the order of the small 
pulsation energy itself. The indication that the trouble is caused by 
advection errors can be glimpsed from Table 1. Here we report results 
obtained with the use of the internal energy equation (Eq. 37) because 
this lets us then monitor total energy conservation. We show both the 
long-term global error e{Etot), defined as the sum of the errors per cycle 
of all the zones divided by E'tot, and the amplitude of the error fluctua- 
tions during the cycle eA{Etot) (also normalized by Etot)- (All the errors 
are below the convergence parameter, 10~^, of the internal iterations of 
the code.) The long term error is indeed seen to be substantially larger 
with the first order advection scheme. When the first order scheme is 
used, the variation of the energy error is dominated by the long term 
trend. 
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Table III. Fourier parameters of the limit cycle solution - magnitude varia- 
tion. (LA: Lagrangean; TE: total energy, IE internal energy Eqs. are solved; 
vL: van Leer, DC: donor cell. 









Scheme 


Po 


Ai 


R21 


R31 


ct>21 


031 


- 


4 


LA 


- 


17.6 


0.41 


0.307 


0.250 


5.815 


5.324 


- 


1 


LA 




IT.G 


0. 12 




{).2.")() 


r).<s:i7 


5.328 


60 


4 


TE 


vL 


17.6 


0.41 


0.299 


0.229 


5.817 


5.202 


60 


4 


IE 


vL 


17.6 


0.41 


0.300 


0.229 


5.825 


5.212 


60 


1 


TE 


vL 


17.6 


0.42 


0.301 


0.248 


5.784 


5.226 


60 


1 


IE 


vL 


17.6 


0.42 


0.301 


0.249 


5.790 


5.227 


40 


4 


TE 


vL 


17.5 


0.40 


0.276 


0.211 


5.878 


5.089 


40 


4 


IE 


vL 


17.5 


0.38 


0.282 


0.209 


5.851 


5.162 


40 


4 


TE 


DC 


18.7 


0.64 


0.364 


0.144 


0.049 


6.082 


40 


4 


IE 


DC 


18.0 


0.31 


0.209 


0.171 


5.712 


4.461 


30 


4 


TE 


vL 


17.3 


0.32 


0.241 


0.151 


5.792 


4.861 


30 


4 


IE 


vL 


17.3 


0.30 


0.248 


0.143 


5.778 


4.990 



The solution to the numerieal instability is to make the motion of 
the inner envelope Lagrangean, as Feuchtinger & Dorfi (1994) have also 
found. The selection of the location of the switching point Jl between 
adaptive and Lagrangean mesh is important, though, as Table 1 shows. 
An optimal switch ji should be as far out in the envelope as possible, 
but above the He II ionization (at T > 5 - 10 x 10^ K). 

Unfortunately no suitable analytical results exist for a stellar pul- 
sator model against which one could check the codes. We therefore 
have to rely on the robustness of the results with respect to numerical 
parameters on the one hand, and to the agreement that we can achieve 
with the observational data. Although it is important that the pul- 
sations should exhibit long term stability (total energy and pulsation 
energies constant), this, by itself, is not a guarantee that the amplitude 
of the cycle is in fact correct. 

In order to compare the robustness of the computed limit cycle 
pulsation we display their Fourier parameters, for the radial velocity of 
the photo-sphere in Table II and for the stellar magnitude in Table III. 

First wc note that the calculations with the new Lagrangean code 
confirm the results of Moskalik, Buchler and Marom (1992) which, we 
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recall, gave Fourier parameters of the bump Cepheid models that were 
rather close to the observed values. 

Second, it immediately stands out that it is imperative to use a 
second order scheme for the fluxes, confirming the results of our test 
examples. We note parenthetically that in a sequence of Cepheid models 
that we calculated the donor cell fluxes produced results with trends 
in the amplitude vs. period behavior that were completely wrong, even 
opposite to those both of the higher order scheme and of the purely La- 
grangean scheme. For the models calculated with Donor cell advection 
scheme, the shifts of amplitudes are even opposite for the two choices 
of handling the energy. 

Third, the Fourier parameters exhibit a reasonable robustness when 
the spatial switchover from Lagrangean to adaptive mesh occurs in the 
range j£=40-60. Values of iL=30 give already unsatisfactory results. 
Note also that when the internal energy (IE) equation is solved the 
Fourier parameters are more sensitive to jl because of the increased ad- 
vection errors. This is especially true for the third order radial velocity 
coefficients R31 and ^31. This again confirms the numerical superiority 
of the total energy over the internal energy equation. 

Fourth, the adaptive code although requiring some pseudo- viscosity, 
lets one reduce pseudo-viscous dissipation below values that are possi- 
ble with the Lagrangean code. The increase of the pulsation amplitude 
with Cq clearly shows that this amplitude is largely determined by 
the artificial viscosity independently of the choice of grid. This again 
illustrates the well known problem that when we reduce the numerical 
dissipation below the 'standard' values, the pulsation amplitudes get 
larger than the observed ones (for a discussion cf. e.g. Kovacs 1990). 
It must be considered a lucky coincidence that, since the pioneering 
work of Christy in the 1960s, the use of the 'standard' vNR pseudo- 
viscosity has yielded pulsation amplitudes in rough agreement with the 
observed ones. In order to make progress it will be necessary to replace 
the pseudo- viscous dissipation with the correct physical one, most likely 
eddy viscosity. It is very likely that the violent shocks induce turbulent 
motions, especially in the relatively 'soft' (small F) partially ionized 
regions. But to even get a good understanding of where in the star 
and when in the pulsation phase such turbulent dissipation plays a 
role, and to get even a crude idea of itsw magnitude, one needs a 3D 
hydrodynamic simulation. 

Finally, we note that the great advantages of the adaptive numerical 
calculations is the smoothness of the results. In Figure 7. the tempera- 
ture variations of the outer 60 shells are shown both for the Lagrangean 
and for the adaptive results (in both cases as a function of the respective 
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grid indices). The resolution of the spatial temperature structure is 
clearly superior by the adaptive mesh. 

In Fig. 8 we present the fluid velocity variations of different layers of 
the model (spread out vertically for clarity), both for the Lagrangean 
and for the adaptive results. The plotted masses have been selected 
at the values MJ^, Mq^, Mqq . . . Mj'^Q of the Langrangean model. In 
order to compare we have had to interpolate the velocities from the 
adaptive code for these given values of the masses. The adaptive code 
is seen to provide significantly smoother curves than the Lagrangean 
The relatively sharp features on some of the velocity curves indicate 
the position of the hydrogen partial ionization front. 

Finally, Figure 9 shows the improvement in the quality of the light- 
curve that is obtained with the adaptive mesh. 



4. Conclusion 

We have described an adaptive hydrocode which has been specifically 
developed for stellar pulsation calculations. The code has been tested 
for shock problems with known analytical solutions. Purposedly, to 
mimic more realistic calculations in which few mesh points are available 
for the shock by itself, a relatively coarse mesh has been used. These 
tests lead to the following conclusions: 

- The first order schemes (donor cell) do not give satisfactory agreement 
between the numerical and analytical solution. At least a second order 
{e.g. van Leer) scheme is necessary. 

- The solution of the total energy equation is preferable to the use of 
an internal energy equation because the numerical deterioration of the 
results due to the advection errors is reduced. 

- A careful application of the adaptive grid technique can provide 
superior results for tough shock problems even with a relatively few 
grid points. 

It is of course much easier and much faster to run a Lagrangean code 
than an adaptive one. In particular, the Stellingwerf Lagrangean relax- 
ation method is very efficient when periodic cycles are to be computed 
as for classical Cepheids or RR Lyrae. The fact that the adaptive mesh 
can be switched on smoothly starting from a Lagrangean description 
is therefore a very useful feature. 

We have applied the code to a typical Cepheid model pulsator. The 
temporal variations of the spatial structure of the temperature are well 
tracked with the adaptive mesh. This results in much smoother velocity 
and luminosity curves than can be obtained with the Lagrangean code. 
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